(library (grid-gen Lax-Friedrichs)

(export Lax-Friedrichs:dF Lax-Friedrichs:dFl Lax-Friedrichs:vm)

(import (guile)
        (grid-gen main)
        (grid-gen math)
        (grid-gen symath)
)

(define set-lr (lambda (v vl vr)
  (lambda ar
    (if (not (null? ar))
        (cond ((eq? (car ar) 'vl) vl)
              ((eq? (car ar) 'vr) vr)
              (#t (apply v ar))
        )
        (apply v ar)
    )
  )
))

(define Lax-Friedrichs:lam (lambda (d)
  (lambda (vm)
    (abs2 ((PHYS:chrs d) vm))
  )
))


(define Lax-Friedrichs:vm (lambda (vl vr)
  (lambda args
    (cond ((eq? (car args) 'vl) vl)
          ((eq? (car args) 'vr) vr)
          (#t `(* 0.5 (+ ,(apply vl args) ,(apply vr args))))
    )
  )
))

(define Lax-Friedrichs:dFl (lambda (d gr i)
  (let* ((vl0 (gr 'var-l i d))
         (vr0 (gr 'var-r i d))
         
         (vl (set-lr vl0 vl0 vr0))
         (vr (set-lr vr0 vl0 vr0))

         (ul (arr-max ((ADV:lam d) vl)))
         (ur (arr-max ((ADV:lam d) vr)))

         (ql (gr 'ql i d))
         (qr (gr 'qr i d))
         
         (amp `(* 0.5 (max ,ul ,ur)))
        )
    (append
       `((,(mul2 amp (add qr (mul2 -1.0 ql)))
           (+ ,(list-ref (list (vr 'u) (vr 'v) (vr 'w)) d) ,amp)
         )
        )
          (map (lamcon `(,(apply math-cvector
                                 (map (lamcon 0.0)
                                      (iota1 (PHYS:qlen (gr 'ndims)))
                                 )
                          )
                          0.0
                        )
               )
               (iota1 (- (PHYS:qlen (gr 'ndims)) 2))
          )
       `((,(mul2 amp (add qr (mul2 -1.0 ql)))
           (- ,(list-ref (list (vl 'u) (vl 'v) (vl 'w)) d) ,amp)
         )
        )
    )
  )
))

(define Lax-Friedrichs:dF (lambda (d gr i)
  (mul2 0.5 (apply add (map car (Lax-Friedrichs:dFl d gr i))))
))


)

(import (grid-gen main))
(set-comp-order 1)
(set! ADV:dF Lax-Friedrichs:dF)
(set! ADV:dFl Lax-Friedrichs:dFl)
(set! ADV:vm Lax-Friedrichs:vm)
(set! ADV:lam Lax-Friedrichs:lam)
(set! ADV:dF1
      (lambda (d gr i)
        (apply math-cvector
               (map (lamcon 0.0) (iota1 (PHYS:qlen (gr 'ndims))))
        )
      )
)
